#include <fstream>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <cmath>
#include "TransWorker.h"
#include "tinyxml2.h"

using namespace std;
using namespace tinyxml2;

// function for edge sort
bool CompareEdge(vector<int> &vEdge1, vector<int> &vEdge2)
{
    if (vEdge1[0] < vEdge2[0]){
        return true;
    } else if (vEdge1[0] == vEdge2[0]){
        if (vEdge1[1] < vEdge2[1]){
            return true;
        } else {
            return false;
        }
    } else {
        return false;
    }
}

bool TransWorker::fnParseCmdArgu(int argc, char *argv[], string &sAedt, string &sMesh)
{
    if (argc < 3) return false;

    bool bAedt = false;
    bool bMesh = false;
    for (int n = 1; n < 3; ++n) {
        auto sFile = string(argv[n]);
        auto i = sFile.length();
        if (sFile.substr(i - 5, 5) == ".aedt") {
            sAedt = sFile;
            bAedt = true;
        } else if (sFile.substr(i - 7, 7) == ".ngmesh") {
            sMesh = sFile;
            bMesh = true;
        }
    }
    if (!bAedt) {
        cout << "Error : Not input aedt file." << endl;
        return false;
    }
    if (!bMesh) {
        cout << "Error : Not input ngmesh file." << endl;
        return false;
    }
    cout << "Aedt file: " << sAedt << endl;
    cout << "Mesh file: " << sMesh << endl;
    return true;
}

bool TransWorker::fnLoadAedtFile(string &sName)
{
    ifstream fp(sName);
    if (!fp.is_open()) return false;

    stringstream buff;
    buff << fp.rdbuf();
    string sFile(buff.str());
    buff.str("");
    fp.close();

    // load material
    auto iHead = sFile.find("$begin 'Materials'");
    auto iTail = sFile.find("$end 'Materials'");
    auto sTxt = sFile.substr(iHead, iTail - iHead);
    fnLoadMaterial(sTxt);

    // load geometry
    iHead = sFile.find("$begin 'ToplevelParts'");
    iTail = sFile.find("$end 'ToplevelParts'");
    sTxt = sFile.substr(iHead, iTail - iHead);
    fnLoadGeometry(sTxt);

    // load boundary
    iHead = sFile.find("$begin 'Boundaries'");
    iTail = sFile.find("$end 'Boundaries'");
    sTxt = sFile.substr(iHead, iTail - iHead);
    fnLoadBoundary(sTxt);

    // load solver setting
    iHead = sFile.find("$begin 'SolveSetups'");
    iTail = sFile.find("$end 'SolveSetups'");
    sTxt = sFile.substr(iHead, iTail - iHead);
    fnLoadSolvPara(sTxt);

    // load port voltage
    iHead = sFile.find("$begin 'Solutions'");
    iTail = sFile.find("$end 'Solutions'");
    sTxt = sFile.substr(iHead, iTail - iHead);
    fnLoadSolution(sTxt);

    // link geometry and material
    bool bOk = fnLinkGeomMedi();
    return bOk;
}

bool TransWorker::fnLoadMeshFile(string &sName)
{
    ifstream fp(sName);
    if (!fp.is_open()) return false;

    if (sName.find("initial.ngmesh") == string::npos) {
        m_para.iAdapMaxNum = 0;
    }
    int iNumBody = 0;
    int iNumFace = 0;
    int iNumNode = 0;
    int iNumTria = 0;
    int iNumTetr = 0;
    string sLine, sKey;
    stringstream iss;
    while (!fp.eof()) {
        getline(fp, sLine);
        if (sLine.find("nbodies") != string::npos) {
            iss.str(sLine);
            iss >> sLine >> iNumBody;
            fp >> sLine >> iNumFace;
            fp >> sLine >> iNumNode;
            fp >> sLine >> iNumTria;
            fp >> sLine >> iNumTetr;
            break;
        }
    }
    if (iNumNode < 1 || iNumTetr < 1) return false;
    // load the body information
    int iBodyID, i, j, m, n;
    while (!fp.eof()) {
        getline(fp, sLine);
        if (sLine.find("begin_body_data") != string::npos) {
            for (i = 0; i < iNumBody; ++i) {
                fp >> sLine >> iBodyID >> sKey >> sKey >> sKey >> n >> sKey >> n >> sKey >> n;
                getline(fp, sLine);
                auto iter = m_mapGeom.find(iBodyID);
                if (iter != m_mapGeom.end()) {
                    fp >> sLine;
                    for (j = 0; j < n; ++j) {
                        fp >> m;
                        m_vGeom[iter->second].vFaceID.push_back(m);
                    }
                }
                for (j = 0; j < 7; ++j) getline(fp, sLine);
            }
            break;
        }
    }
    // map the boundary and face
    map<int, int> mapFace;
    fnSetBoundFace(mapFace);

    m_vR.resize(iNumNode);
    m_vTria.reserve(iNumTria);
    m_vTetr.reserve(iNumTetr);
    // load the points
    while (!fp.eof()) {
        getline(fp, sLine);
        if (sLine.find("begin_point_data") != string::npos) {
            int n;
            string sKey;
            for (int i = 0; i < iNumNode; ++i) {
                fp >> sKey >> n >> sLine >> m_vR[i].x >> m_vR[i].y >> m_vR[i].z;
            }
            break;
        }
    }
    // load the facet element
    m_vTria.clear();
    size_t indx;
    while (!fp.eof()) {
        getline(fp, sLine);
        if (sLine.find("begin_facet_data") != string::npos) {
            string sKey;
            Tria tri{};
            for (int i = 0; i < iNumFace; ++i) {
                fp >> sLine >> m >> sKey >> m >> sKey >> n;
                getline(fp, sLine);
                indx = sLine.find("face_ids") + 9;
                iss.clear();
                iss.str(sLine.substr(indx));
                tri.iBoud = -1;
                while (!iss.eof() && iss >> m) {
                    auto iter = mapFace.find(m);
                    if (iter == mapFace.end()) continue;
                    if (-1 == tri.iBoud) {
                        tri.iBoud = iter->second;
                    } else {
                        if (m_vBoud[tri.iBoud].sType != WAVE_PORT &&
                            m_vBoud[tri.iBoud].sType != LUMP_PORT) {
                            tri.iBoud = iter->second;
                        }
                    }
                }
                
                for (m = 0; m < 6; ++m ) getline(fp, sLine);
                for (int j = 0; j < n; ++j) {
                    fp >> sLine >> m >> sKey >> m >> sKey >> m;
                    fp >> sKey >> tri.node[0] >> tri.node[1] >> tri.node[2];
                    if (tri.iBoud != -1) m_vTria.push_back(tri);
                }
            }
            break;
        }
    }
    // load the volume element
    m_vTetr.clear();
    while (!fp.eof()) {
        getline(fp, sLine);
        if (sLine.find("begin_vol_element_data") != string::npos) {
            Tetr tet{};
            for (int i = 0; i < iNumTetr; ++i) {
                fp >> sKey >> n >> sLine >> tet.iBody >> sKey >> m >> sLine;
                fp >> tet.node[0] >> tet.node[1] >> tet.node[2] >> tet.node[3] >> sKey >> n;
                m_vTetr.push_back(tet);
            }
            break;
        }
    }
    fp.close();
    // set Tetrahedra material
    fnTetrMaterial();
    // get the triangle on each boundary
    fnGetBoundTria();
    // link triangle and Tetrahedra
    fnLinkTriaTetr();

    return true;
}

bool TransWorker::fnGetPortCoord()
{
    Node r1, r2;
    double dLen;

    vector<int> vWire;
    for (auto &lump : m_vLump) {
        fnGetWireFrame(lump, vWire);
        if (vWire.size() == 4) {
            lump.sType = "Rectangle";
            r1 = m_vR[vWire[1]] - m_vR[vWire[0]];
            r2 = m_vR[vWire[3]] - m_vR[vWire[0]];
            if (r1.z < 0) r1 = -r1;
            if (r2.z < 0) r2 = -r2;
            if (r1.z > r2.z) {
                lump.ei = r1;
                lump.a = get_length(r2);
                lump.b = get_length(r1);
            } else {
                lump.ei = r2;
                lump.a = get_length(r1);
                lump.b = get_length(r2);
            }
            lump.ei /= lump.b;
        } else {
            lump.sType = "Coaxial";
            lump.a = 1.0e20;
            lump.b = 0.0;
            for (auto &wire : vWire) {
                r1 = m_vR[wire] - lump.rc;
                dLen = get_length(r1);
                if (lump.a > dLen) lump.a = dLen;
                if (lump.b < dLen) lump.b = dLen;
            }
        }

        if (abs(lump.ei.x) < 1.0e-15) lump.ei.x = 0;
        if (abs(lump.ei.y) < 1.0e-15) lump.ei.y = 0;
        if (abs(lump.ei.z) < 1.0e-15) lump.ei.z = 0;
    }
    for (auto &port : m_vPort) {
        fnGetWireFrame(port, vWire);
        if (vWire.size() == 4) {
            port.sType = "Rectangle";
            r1 = m_vR[vWire[1]] - m_vR[vWire[0]];
            r2 = m_vR[vWire[3]] - m_vR[vWire[0]];
            if (r1.z < 0) r1 = -r1;
            if (r2.z < 0) r2 = -r2;
            if (r1.z > r2.z) {
                port.ei = r1;
                port.a = get_length(r2);
                port.b = get_length(r1);
            } else {
                port.ei = r2;
                port.a = get_length(r1);
                port.b = get_length(r2);
            }
            port.ei /= port.b;
        } else {
            port.sType = "General";
            port.a = 1.0e20;
            port.b = 0.0;
            for (auto &wire : vWire) {
                r1 = m_vR[wire] - port.rc;
                dLen = get_length(r1);
                if (port.a > dLen) port.a = dLen;
                if (port.b < dLen) port.b = dLen;
            }
        }
        if (abs(port.ei.x) < 1.0e-15) port.ei.x = 0;
        if (abs(port.ei.y) < 1.0e-15) port.ei.y = 0;
        if (abs(port.ei.z) < 1.0e-15) port.ei.z = 0;
    }
    return true;
}

void TransWorker::fnLoadMaterial(string &sTxt)
{
    stringstream iss(sTxt);
    string sLine;

    size_t iLen, indx;
    string sKey;
    m_vMedi.clear();
    m_vMedi.emplace_back();  // set index 0 as PEC
    m_vMedi[0].sName = "PEC";
    m_vMedi[0].dSgm = 1.0e20;
    getline(iss, sLine);
    Material medi;
    while (!iss.eof()) {
        getline(iss, sLine);
        if (sLine.find("$begin") == string::npos) break;
        indx = sLine.find("'")+1;
        iLen = sLine.length()-indx-1;
        medi.sName = sLine.substr(indx, iLen);
        medi.dEps = 1.0;
        medi.dTan = 0.0;
        medi.dSgm = 0.0;
        medi.dMiu = 1.0;
        medi.dTanM = 0.0;
        medi.iTensor = 0;
        while (true) {
            getline(iss, sLine);
            if (sLine.find("permittivity") != string::npos) {
                if (medi.iTensor == 0 && sLine.find("$begin") == string::npos) {
                    iLen = sLine.length() - 1;
                    indx = sLine.find('=') + 2;
                    sKey = sLine.substr(indx, iLen - indx);
                    medi.dEps = stod(sKey);
                } else {
                    medi.iTensor = 1;
                    vector<double> vEps;
                    while (true) {
                        getline(iss, sLine);
                        if (sLine.find("component") != string::npos) {
                            iLen = sLine.length() - 1;
                            indx = sLine.find('=') + 2;
                            sKey = sLine.substr(indx, iLen - indx);
                            vEps.emplace_back(stod(sKey));

                        } else if(sLine.find("$end") != string::npos) {
                            break;
                        }
                    }
                    if (vEps.size() == 3) {
                        medi.vEpsRe[0][0] = vEps[0];
                        medi.vEpsRe[1][1] = vEps[1];
                        medi.vEpsRe[2][2] = vEps[2];
                    } else if (vEps.size() == 9) {
                        medi.vEpsRe[0][0] = vEps[0];
                        medi.vEpsRe[0][1] = vEps[1];
                        medi.vEpsRe[0][2] = vEps[2];
                        medi.vEpsRe[1][0] = vEps[3];
                        medi.vEpsRe[1][1] = vEps[4];
                        medi.vEpsRe[1][2] = vEps[5];
                        medi.vEpsRe[2][0] = vEps[6];
                        medi.vEpsRe[2][1] = vEps[7];
                        medi.vEpsRe[2][2] = vEps[8];
                    }
                }
            } else if (sLine.find("dielectric_loss_tangent") != string::npos) {
                if (medi.iTensor == 0 && sLine.find("$begin") == string::npos) {
                    iLen = sLine.length() - 1;
                    indx = sLine.find('=') + 2;
                    sKey = sLine.substr(indx, iLen - indx);
                    medi.dTan = stod(sKey);
                } else {
                    vector<double> vEps;
                    while (true) {
                        getline(iss, sLine);
                        if (sLine.find("component") != string::npos) {
                            iLen = sLine.length() - 1;
                            indx = sLine.find('=') + 2;
                            sKey = sLine.substr(indx, iLen - indx);
                            vEps.emplace_back(stod(sKey));

                        } else if (sLine.find("$end") != string::npos) {
                            break;
                        }
                    }
                    if (vEps.size() == 3) {
                        medi.vEpsIm[0][0] = vEps[0];
                        medi.vEpsIm[1][1] = vEps[1];
                        medi.vEpsIm[2][2] = vEps[2];
                    } else if (vEps.size() == 9) {
                        medi.vEpsIm[0][0] = vEps[0];
                        medi.vEpsIm[0][1] = vEps[1];
                        medi.vEpsIm[0][2] = vEps[2];
                        medi.vEpsIm[1][0] = vEps[3];
                        medi.vEpsIm[1][1] = vEps[4];
                        medi.vEpsIm[1][2] = vEps[5];
                        medi.vEpsIm[2][0] = vEps[6];
                        medi.vEpsIm[2][1] = vEps[7];
                        medi.vEpsIm[2][2] = vEps[8];
                    }
                }
            } else if (sLine.find("permeability") != string::npos) {
                if (medi.iTensor == 0 && sLine.find("$begin") == string::npos) {
                    iLen = sLine.length() - 1;
                    indx = sLine.find('=') + 2;
                    sKey = sLine.substr(indx, iLen - indx);
                    medi.dMiu = stod(sKey);
                } else {
                    vector<double> vMiu;
                    while (true) {
                        getline(iss, sLine);
                        if (sLine.find("component") != string::npos) {
                            iLen = sLine.length() - 1;
                            indx = sLine.find('=') + 2;
                            sKey = sLine.substr(indx, iLen - indx);
                            vMiu.emplace_back(stod(sKey));

                        } else if (sLine.find("$end") != string::npos) {
                            break;
                        }
                    }
                    if (vMiu.size() == 3) {
                        medi.vMiuRe[0][0] = vMiu[0];
                        medi.vMiuRe[1][1] = vMiu[1];
                        medi.vMiuRe[2][2] = vMiu[2];
                    } else if (vMiu.size() == 9) {
                        medi.vMiuRe[0][0] = vMiu[0];
                        medi.vMiuRe[0][1] = vMiu[1];
                        medi.vMiuRe[0][2] = vMiu[2];
                        medi.vMiuRe[1][0] = vMiu[3];
                        medi.vMiuRe[1][1] = vMiu[4];
                        medi.vMiuRe[1][2] = vMiu[5];
                        medi.vMiuRe[2][0] = vMiu[6];
                        medi.vMiuRe[2][1] = vMiu[7];
                        medi.vMiuRe[2][2] = vMiu[8];
                    }
                }
            } else if (sLine.find("magnetic_loss_tangent") != string::npos) {
                if (medi.iTensor == 0 && sLine.find("$begin") == string::npos) {
                    iLen = sLine.length() - 1;
                    indx = sLine.find('=') + 2;
                    sKey = sLine.substr(indx, iLen - indx);
                    medi.dTanM = stod(sKey);
                } else {
                    vector<double> vMiu;
                    while (true) {
                        getline(iss, sLine);
                        if (sLine.find("component") != string::npos) {
                            iLen = sLine.length() - 1;
                            indx = sLine.find('=') + 2;
                            sKey = sLine.substr(indx, iLen - indx);
                            vMiu.emplace_back(stod(sKey));

                        } else if (sLine.find("$end") != string::npos) {
                            break;
                        }
                    }
                    if (vMiu.size() == 3) {
                        medi.vMiuIm[0][0] = vMiu[0];
                        medi.vMiuIm[1][1] = vMiu[1];
                        medi.vMiuIm[2][2] = vMiu[2];
                    } else if (vMiu.size() == 9) {
                        medi.vMiuIm[0][0] = vMiu[0];
                        medi.vMiuIm[0][1] = vMiu[1];
                        medi.vMiuIm[0][2] = vMiu[2];
                        medi.vMiuIm[1][0] = vMiu[3];
                        medi.vMiuIm[1][1] = vMiu[4];
                        medi.vMiuIm[1][2] = vMiu[5];
                        medi.vMiuIm[2][0] = vMiu[6];
                        medi.vMiuIm[2][1] = vMiu[7];
                        medi.vMiuIm[2][2] = vMiu[8];
                    }
                }
            } else if (sLine.find("_conductivity=") != string::npos) {
                ;   // do nothing, skip "thermal_conductivity"
            } else if (sLine.find("conductivity=") != string::npos) {
                iLen = sLine.length() - 1;
                indx = sLine.find('=') + 2;
                sKey = sLine.substr(indx, iLen - indx);
                medi.dSgm = stod(sKey);
            } else if (sLine.find(medi.sName) != string::npos) {
                break;
            }
        }
        if(medi.sName == "perfect conductor" || medi.sName == "pec") continue;
        if (medi.iTensor == 1) {
            double temp[3][3];
            for (int i = 0; i < 3; ++i) {
                for (int j = 0; j < 3; ++j) {
                    temp[i][j] = medi.vEpsRe[i][0] * medi.vEpsIm[0][j] + medi.vEpsRe[i][1] * medi.vEpsIm[1][j] + medi.vEpsRe[i][2] * medi.vEpsIm[2][j];
                }
            }
            for (int i = 0; i < 3; ++i) {
                for (int j = 0; j < 3; ++j) {
                    medi.vEpsIm[i][j] = -temp[i][j];
                }
            }
            for (int i = 0; i < 3; ++i) {
                for (int j = 0; j < 3; ++j) {
                    temp[i][j] = medi.vMiuRe[i][0] * medi.vMiuIm[0][j] + medi.vMiuRe[i][1] * medi.vMiuIm[1][j] + medi.vMiuRe[i][2] * medi.vMiuIm[2][j];
                }
            }
            for (int i = 0; i < 3; ++i) {
                for (int j = 0; j < 3; ++j) {
                    medi.vMiuIm[i][j] = -temp[i][j];
                }
            }
            int first_ = medi.sName.find_first_of("_");
            int last_ = medi.sName.find_last_of("_");
            auto pmlType = medi.sName.substr(first_ + 1, last_ - first_ - 1);
            if (pmlType == "XZ") {
                Material mediXY;
                mediXY.iTensor = 1;
                mediXY.sName = medi.sName.substr(0, first_ + 1) + "XY" + medi.sName.substr(last_, medi.sName.size() - last_);
                for (int i = 0; i < 3; ++i) {
                    for (int j = 0; j < 3; ++j) {
                        mediXY.vEpsRe[i][j] = medi.vEpsRe[i][j];
                        mediXY.vEpsIm[i][j] = medi.vEpsIm[i][j];
                        mediXY.vMiuRe[i][j] = medi.vMiuRe[i][j];
                        mediXY.vMiuIm[i][j] = medi.vMiuIm[i][j];
                        if (i == 1 && j == 1) {
                            mediXY.vEpsRe[i][j] = medi.vEpsRe[i + 1][j + 1];
                            mediXY.vEpsIm[i][j] = medi.vEpsIm[i + 1][j + 1];
                            mediXY.vMiuRe[i][j] = medi.vMiuRe[i + 1][j + 1];
                            mediXY.vMiuIm[i][j] = medi.vMiuIm[i + 1][j + 1];
                        } else if (i == 2 && j == 2) {
                            mediXY.vEpsRe[i][j] = medi.vEpsRe[i - 1][j - 1];
                            mediXY.vEpsIm[i][j] = medi.vEpsIm[i - 1][j - 1];
                            mediXY.vMiuRe[i][j] = medi.vMiuRe[i - 1][j - 1];
                            mediXY.vMiuIm[i][j] = medi.vMiuIm[i - 1][j - 1];
                        }
                    }
                }
                m_vMedi.push_back(mediXY);
            } else if (pmlType == "Z") {
                Material mediX;
                mediX.iTensor = 1;
                mediX.sName = medi.sName.substr(0, first_ + 1) + "X" + medi.sName.substr(last_, medi.sName.size() - last_);
                for (int i = 0; i < 3; ++i) {
                    for (int j = 0; j < 3; ++j) {
                        mediX.vEpsRe[i][j] = medi.vEpsRe[i][j];
                        mediX.vEpsIm[i][j] = medi.vEpsIm[i][j];
                        mediX.vMiuRe[i][j] = medi.vMiuRe[i][j];
                        mediX.vMiuIm[i][j] = medi.vMiuIm[i][j];
                        if (i == 0 && j == 0) {
                            mediX.vEpsRe[i][j] = medi.vEpsRe[i + 2][j + 2];
                            mediX.vEpsIm[i][j] = medi.vEpsIm[i + 2][j + 2];
                            mediX.vMiuRe[i][j] = medi.vMiuRe[i + 2][j + 2];
                            mediX.vMiuIm[i][j] = medi.vMiuIm[i + 2][j + 2];
                        } else if (i == 2 && j == 2) {
                            mediX.vEpsRe[i][j] = medi.vEpsRe[i - 2][j - 2];
                            mediX.vEpsIm[i][j] = medi.vEpsIm[i - 2][j - 2];
                            mediX.vMiuRe[i][j] = medi.vMiuRe[i - 2][j - 2];
                            mediX.vMiuIm[i][j] = medi.vMiuIm[i - 2][j - 2];
                        }
                    }
                }
                m_vMedi.push_back(mediX);
                Material mediY;
                mediY.iTensor = 1;
                mediY.sName = medi.sName.substr(0, first_ + 1) + "Y" + medi.sName.substr(last_, medi.sName.size() - last_);
                for (int i = 0; i < 3; ++i) {
                    for (int j = 0; j < 3; ++j) {
                        mediY.vEpsRe[i][j] = medi.vEpsRe[i][j];
                        mediY.vEpsIm[i][j] = medi.vEpsIm[i][j];
                        mediY.vMiuRe[i][j] = medi.vMiuRe[i][j];
                        mediY.vMiuIm[i][j] = medi.vMiuIm[i][j];
                        if (i == 1 && j == 1) {
                            mediY.vEpsRe[i][j] = medi.vEpsRe[i + 1][j + 1];
                            mediY.vEpsIm[i][j] = medi.vEpsIm[i + 1][j + 1];
                            mediY.vMiuRe[i][j] = medi.vMiuRe[i + 1][j + 1];
                            mediY.vMiuIm[i][j] = medi.vMiuIm[i + 1][j + 1];
                        } else if (i == 2 && j == 2) {
                            mediY.vEpsRe[i][j] = medi.vEpsRe[i - 1][j - 1];
                            mediY.vEpsIm[i][j] = medi.vEpsIm[i - 1][j - 1];
                            mediY.vMiuRe[i][j] = medi.vMiuRe[i - 1][j - 1];
                            mediY.vMiuIm[i][j] = medi.vMiuIm[i - 1][j - 1];
                        }
                    }
                }
                m_vMedi.push_back(mediY);
            }
        }
        m_vMedi.push_back(medi);
    }
    // move vacuum to the second
    for (int i = 2; i < m_vMedi.size(); ++i) {
        if (m_vMedi[i].sName == "vacuum") {
            medi = m_vMedi[i];
            m_vMedi[i] = m_vMedi[1];
            m_vMedi[1] = medi;
        }
    }
}

void TransWorker::fnLoadGeometry(string &sTxt)
{
    stringstream iss(sTxt);
    string sLine, sKey;

    size_t iLen, indx;
    m_vGeom.clear();
    while (!iss.eof()) {
        getline(iss, sLine);
        if (sLine.find("$begin 'GeometryPart'") != string::npos) {
            Geometry geom;
            while (sLine.find("$end 'Attributes'") == string::npos) {
                getline(iss, sLine);
                if (sLine.find("Name=") != string::npos) {
                    indx = sLine.find('=') + 2;
                    iLen = sLine.length() - indx - 1;
                    geom.sName = sLine.substr(indx, iLen);
                } else if (sLine.find("SurfaceMaterialValue") != string::npos) {
                    ;   // do nothing, this sentence is very important
                } else if (sLine.find("MaterialValue=") != string::npos) {
                    indx = sLine.find('=') + 3;
                    iLen = sLine.length() - indx - 2;
                    geom.sMedi = sLine.substr(indx, iLen);
                    sKey.resize(geom.sMedi.length());
                    transform(geom.sMedi.begin(), geom.sMedi.end(), sKey.begin(), ::tolower);
                    if(geom.sMedi.size() == 0) geom.sMedi = "vacuum";
                    if (sKey == "perfect conductor" || sKey == "pec") geom.sMedi = "PEC";
                } else if (sLine.find("SolveInside=") != string::npos) {
                    indx = sLine.find('=') + 1;
                    sKey = sLine.substr(indx);
                    if (sKey == "true") {
                        geom.bSolvInside = true;
                    } else {
                        geom.bSolvInside = false;
                    }
                }
            }
            getline(iss, sLine);
            while (sLine.find("$end 'Operation'") == string::npos) {
                getline(iss, sLine);
                if (sLine.find("ParentPartID") != string::npos) {
                    indx = sLine.find('=') + 1;
                    sKey = sLine.substr(indx);
                    geom.iObjID = stoi(sKey);
                }
            }
            while (sLine.find("$end 'GeometryPart'") == string::npos) {
                getline(iss, sLine);
            }
            m_vGeom.push_back(geom);
        }
    }
    for (int i = 0; i < m_vGeom.size(); ++i) {
        if (m_vGeom[i].sName.find("PML") != string::npos) {
            int firstName = m_vGeom[i].sName.find_first_of("_");
            int lastName = m_vGeom[i].sName.find_last_of("_");
            string lastNum = m_vGeom[i].sName.substr(lastName + 1);
            int firstMedia = m_vGeom[i].sMedi.find_first_of("_");
            int lastMedia = m_vGeom[i].sMedi.find_last_of("_");
            string pmlType = m_vGeom[i].sMedi.substr(firstMedia + 1, lastMedia - firstMedia - 1);
            if (pmlType == string("Z") && (lastNum == string("4") || lastNum == string("6"))) {
                m_vGeom[i].sMedi = m_vGeom[i].sMedi.substr(0, firstMedia+1) + "X" + m_vGeom[i].sMedi.substr(lastMedia, m_vGeom[i].sMedi.size() - lastMedia);
            } else if (pmlType == string("Z") && (lastNum == string("3") || lastNum == string("5"))) {
                m_vGeom[i].sMedi = m_vGeom[i].sMedi.substr(0, firstMedia + 1) + "Y" + m_vGeom[i].sMedi.substr(lastMedia, m_vGeom[i].sMedi.size() - lastMedia);
            } else if (pmlType == string("XZ") && (lastNum == string("7") || lastNum == string("9") || lastNum == string("13"))) {
                m_vGeom[i].sMedi = m_vGeom[i].sMedi.substr(0, firstMedia + 1) + "YZ" + m_vGeom[i].sMedi.substr(lastMedia, m_vGeom[i].sMedi.size() - lastMedia);
            } else if (pmlType == string("YZ") && (lastNum == string("10") || lastNum == string("12") || lastNum == string("14"))) {
                m_vGeom[i].sMedi = m_vGeom[i].sMedi.substr(0, firstMedia + 1) + "XZ" + m_vGeom[i].sMedi.substr(lastMedia, m_vGeom[i].sMedi.size() - lastMedia);
            } else if (pmlType == string("YZ") && (lastNum == string("15") || lastNum == string("16") || lastNum == string("17") || lastNum == string("18"))) {
                m_vGeom[i].sMedi = m_vGeom[i].sMedi.substr(0, firstMedia + 1) + "XY" + m_vGeom[i].sMedi.substr(lastMedia, m_vGeom[i].sMedi.size() - lastMedia);
            } else if (pmlType == string("YZ") && (lastNum == string("11"))) {
                m_vGeom[i].sMedi = m_vGeom[i].sMedi.substr(0, firstMedia + 1) + "YZ" + m_vGeom[i].sMedi.substr(lastMedia, m_vGeom[i].sMedi.size() - lastMedia);
            } else if (pmlType == string("XZ") && (lastNum == string("8"))) {
                m_vGeom[i].sMedi = m_vGeom[i].sMedi.substr(0, firstMedia + 1) + "XZ" + m_vGeom[i].sMedi.substr(lastMedia, m_vGeom[i].sMedi.size() - lastMedia);
            }
        }
    }
    // create the map from geometry ID to indx
    m_mapGeom.clear();
    for (int i = 0; i < m_vGeom.size(); ++i) {
        m_mapGeom.insert(pair<int, int>(m_vGeom[i].iObjID, i));
    }
}

void TransWorker::fnLoadBoundary(string &sTxt)
{
    stringstream iss(sTxt);
    string sLine, sEnd, sKey;
    m_vBoud.clear();

    size_t iLen, indx;
    getline(iss, sLine);
    getline(iss, sLine);
    while (sLine.find("$begin") == string::npos) {
        getline(iss, sLine);
    }
    while (!iss.eof()) {
        Boundary boud;
        indx = sLine.find('\'');
        boud.sName = sLine.substr(indx);
        getline(iss, sLine);
        indx = sLine.find('=') + 1;
        sKey = sLine.substr(indx);
        boud.id = stoi(sKey);
        getline(iss, sLine);
        indx = sLine.find('=') + 1;
        boud.sType = sLine.substr(indx);
        iLen = boud.sType.length() - 2;
        boud.sType = boud.sType.substr(1, iLen);
        sEnd = "$end " + boud.sName;
        while (sLine.find(sEnd) == string::npos) {
            getline(iss, sLine);
            if (sLine.find("Objects") != string::npos) {
                indx = sLine.find('(') + 1;
                iLen = sLine.length() - indx - 1;
                sKey = sLine.substr(indx, iLen);
                indx = sKey.find(',');
                while (indx != string::npos) {
                    boud.vVoluID.push_back(stoi(sKey.substr(0, indx)));
                    sKey = sKey.substr(indx + 1);
                    indx = sKey.find(',');
                }
                boud.vVoluID.push_back(stoi(sKey));
            } else if (sLine.find("Faces") != string::npos) {
                indx = sLine.find('(') + 1;
                iLen = sLine.length() - indx - 1;
                sKey = sLine.substr(indx, iLen);
                indx = sKey.find(',');
                while (indx != string::npos) {
                    boud.vFaceID.push_back(stoi(sKey.substr(0, indx)));
                    sKey = sKey.substr(indx + 1);
                    indx = sKey.find(',');
                }
                boud.vFaceID.push_back(stoi(sKey));
            } else if (sLine.find("$begin 'GeometryPosition'") != string::npos) {
                getline(iss, sLine);
                while (sLine.find("EntityID=") == string::npos) {
                    getline(iss, sLine);
                }
                indx = sLine.find('=') + 1;
                sKey = sLine.substr(indx);
                boud.iEdgeDn = stoi(sKey);
                getline(iss, sLine);
                while (sLine.find("EntityID=") == string::npos) {
                    getline(iss, sLine);
                }
                indx = sLine.find('=') + 1;
                sKey = sLine.substr(indx);
                boud.iEdgeUp = stoi(sKey);
            } else if (sLine.find("Impedance=") != string::npos) {
                indx = sLine.find('=') + 2;
                iLen = sLine.length() - indx - 4;
                sKey = sLine.substr(indx, iLen);
                boud.dZ0 = stod(sKey);
            }
        }
        iLen = boud.sName.length() - 2;
        boud.sName = boud.sName.substr(1, iLen);
        if (boud.sType != "Terminal") m_vBoud.push_back(boud);
        getline(iss, sLine);
    }
    // set boundary type
    for (auto &boud : m_vBoud) {
        if (boud.sType == "Perfect E") {
            boud.sType = "PEC";
        } else if (boud.sType == "Perfect H") {
            boud.sType = "PMC";
        } else if (boud.sType == "Radiation") {
            boud.sType = "ABC";
        } else if (boud.sType == "Lumped Port") {
            boud.sType = LUMP_PORT;
        } else if (boud.sType == "Wave Port") {
            boud.sType = WAVE_PORT;
        }
    }
}

void TransWorker::fnLoadSolvPara(string &sTxt)
{
    stringstream iss(sTxt);
    string sLine, sEnd, sKey;

    size_t iLen, indx;
    getline(iss, sLine);
    getline(iss, sLine);
    while (!iss.eof()) {
        indx = sLine.find("$begin");
        if (indx != string::npos) break;
        getline(iss, sLine);
    }
    sEnd = "$end" + sLine.substr(indx + 6);
    while (!iss.eof()) {
        if (sLine.find(sEnd) != string::npos) break;
        getline(iss, sLine);
        if (sLine.find("Frequency=") != string::npos) {
            indx = sLine.find('=') + 2;
            iLen = sLine.length() - indx - 3;
            sKey = sLine.substr(indx, iLen);
            indx = sKey.length() - 1;
            if (sKey[indx] == 'K') {
                m_para.sFreqUnit = "KHz";
                sKey = sKey.substr(0, indx);
            } else if (sKey[indx] == 'M') {
                m_para.sFreqUnit = "MHz";
                sKey = sKey.substr(0, indx);
            } else if (sKey[indx] == 'G') {
                m_para.sFreqUnit = "GHz";
                sKey = sKey.substr(0, indx);
            } else {
                m_para.sFreqUnit = "Hz";
            }
            m_para.dFreq0 = stod(sKey);
            m_para.sFreqUnit = sLine.substr(sLine.length() - 4, 3);
        } else if (sLine.find("MaxDeltaS=") != string::npos) {
            indx = sLine.find('=') + 1;
            sKey = sLine.substr(indx);
            m_para.dAdaptDelta = stod(sKey);
        } else if (sLine.find("MaximumPasses=") != string::npos) {
            indx = sLine.find('=') + 1;
            sKey = sLine.substr(indx);
            m_para.iAdapMaxNum = stoi(sKey);
        } else if (sLine.find("MinimumPasses=") != string::npos) {
            indx = sLine.find('=') + 1;
            sKey = sLine.substr(indx);
            m_para.iAdapMinNum = stoi(sKey);
        } else if (sLine.find("MinimumConvergedPasses=") != string::npos) {
            indx = sLine.find('=') + 1;
            sKey = sLine.substr(indx);
            m_para.iAdapPassed = stoi(sKey);
        } else if (sLine.find("PercentRefinement=") != string::npos) {
            indx = sLine.find('=') + 1;
            sKey = sLine.substr(indx);
            m_para.dAdaptRatio = stod(sKey);
        } else if (sLine.find("BasisOrder=") != string::npos) {
            indx = sLine.find('=') + 1;
            sKey = sLine.substr(indx);
            m_para.iOrder = stoi(sKey) + 1;
            if (m_para.iOrder > 2) m_para.iOrder = 2;
        } else if (sLine.find("$begin 'Sweeps'") != string::npos) {
            while (sLine.find("$end 'Sweeps'") == string::npos) {
                getline(iss, sLine);
                if (sLine.find("RangeStart=") != string::npos) {
                    indx = sLine.find('=') + 2;
                    iLen = sLine.length() - indx - 3;
                    sKey = sLine.substr(indx, iLen);
                    indx = sKey.length() - 1;
                    if (sKey[indx] == 'K' || sKey[indx] == 'M' || sKey[indx] == 'G') {
                        sKey = sKey.substr(0, indx);
                    }
                    m_para.dMinFreq = stod(sKey);
                } else if (sLine.find("RangeEnd=") != string::npos) {
                    indx = sLine.find('=') + 2;
                    iLen = sLine.length() - indx - 3;
                    sKey = sLine.substr(indx, iLen);
                    indx = sKey.length() - 1;
                    if (sKey[indx] == 'K' || sKey[indx] == 'M' || sKey[indx] == 'G') {
                        sKey = sKey.substr(0, indx);
                    }
                    m_para.dMaxFreq = stod(sKey);
                } else if (sLine.find("RangeCount=") != string::npos) {
                    indx = sLine.find('=') + 1;
                    sKey = sLine.substr(indx);
                    m_para.iNumFreq = stoi(sKey);
                } else if (sLine.find("RangeType=") != string::npos) {
                    ;  // skip RangeType
                } else if (sLine.find("Type=") != string::npos) {
                    indx = sLine.find('=') + 1;
                    sKey = sLine.substr(indx);
                    if (sKey == "'Interpolating'") {
                        m_para.sSwpType = "Interpolate";
                    } else {
                        m_para.sSwpType = "Discrete";
                    }
                }
            }
        }
    }
}

void TransWorker::fnLoadSolution(string &sTxt)
{
    m_vLump.clear();
    m_vPort.clear();
    for (auto &boud : m_vBoud) {
        if (boud.sType == LUMP_PORT) {
            m_vLump.emplace_back();
            m_vLump.back().id = boud.id;
            m_vLump.back().sName = "lump"+boud.sName;
            m_vLump.back().iEdgeDn = boud.iEdgeDn;
            m_vLump.back().iEdgeUp = boud.iEdgeUp;
            m_vLump.back().dZ0 = boud.dZ0;
        } else if (boud.sType == WAVE_PORT) {
            m_vPort.emplace_back();
            m_vPort.back().id = boud.id;
            m_vPort.back().sName = "port"+boud.sName;
            m_vPort.back().iEdgeDn = boud.iEdgeDn;
            m_vPort.back().iEdgeUp = boud.iEdgeUp;
            m_vPort.back().dZ0 = boud.dZ0;
        }
    }
    map<int, int> mapLump;
    for (int i = 0; i < m_vLump.size(); ++i) {
        mapLump.insert(pair<int, int>(m_vLump[i].id, i));
    }
    map<int, int> mapPort;
    for (int i = 0; i < m_vPort.size(); ++i) {
        mapPort.insert(pair<int, int>(m_vPort[i].id, i));
    }

    stringstream iss(sTxt);
    string sLine, sKey;
    size_t iLen, indx;
    int id;
    PortInfo *port = nullptr;
    while (!iss.eof()) {
        getline(iss, sLine);
        if (sLine.find("SourceEntry") != string::npos) {
            indx = sLine.find("ID=") + 3;
            iLen = sLine.find(",") - indx;
            sKey = sLine.substr(indx, iLen);
            id = stoi(sKey);
            auto iter = mapPort.find(id);
            if (iter != mapPort.end()) {
                port = &(m_vPort[iter->second]);
            } else {
                iter = mapLump.find(id);
                port = &(m_vLump[iter->second]);
            }
            if (sLine.find("Magnitude=") == string::npos) continue;
            indx = sLine.find("Magnitude=") + 11;
            sLine = sLine.substr(indx);
            indx = sLine.find('\'') - 2;
            if (sLine[indx] < 48 || sLine[indx] > 57) {
                sKey = sLine.substr(0, indx);
                port->dMag = stod(sKey);
                sKey = sLine.substr(indx, 2);
                if (sKey == "mW") {
                    port->dMag *= 1.0e-3;
                } else if (sKey == "uW") {
                    port->dMag *= 1.0e-6;
                }
            } else {
                sKey = sLine.substr(0, indx+1);
                port->dMag = stod(sKey);
            }
            port->dMag = sqrt(port->dMag);
            indx = sLine.find("Phase=") + 7;
            sLine = sLine.substr(indx);
            indx = sLine.find('\'') - 3;
            sKey = sLine.substr(0, indx);
            port->dPha = stod(sKey);
            if (sLine.find("rad") != string::npos) {
                 port->dPha *= (180/3.141592653);
            }
        }
    }
}

void TransWorker::fnLoadPostPara(string &sTxt)
{

}

bool TransWorker::fnLinkGeomMedi()
{
    // create the map of material name and index
    map<string, int> mapMedi;
    for (auto i = 0; i < m_vMedi.size(); ++i) {
        mapMedi.insert(pair<string, int>(m_vMedi[i].sName, i));
    }
    // set the link
    for (auto &geo : m_vGeom) {
        auto iter = mapMedi.find(geo.sMedi);
        if (iter == mapMedi.end()) {
            cout << "Error: Gemo. " << geo.sName << " no material." << endl;
            return false;
        } else {
            geo.iMedi = iter->second;
        }
    }
    return true;
}

void TransWorker::fnSetBoundFace(map<int, int> &mapFace)
{
    mapFace.clear();
    for (int i = 0; i < m_vBoud.size(); ++i) {
        // insert the face on boundary
        for (auto &j : m_vBoud[i].vFaceID) {
            mapFace.insert(pair<int, int>(j, i));
        }
        // insert all face on object
        for (auto &j : m_vBoud[i].vVoluID) {
            auto iter = m_mapGeom.find(j);
            if (iter == m_mapGeom.end()) continue;
            for (auto &n : m_vGeom[iter->second].vFaceID) {
                mapFace.insert(pair<int, int>(n, i));
            }
        }
    }
}

void TransWorker::fnGetBoundTria()
{
    auto tmpTri = m_vTria;
    m_vTria.clear();

    for (int i = 0; i < m_vBoud.size(); ++i) {
        m_vBoud[i].iTriaHead = int(m_vTria.size());
        for (auto &tri : tmpTri) {
            if (tri.iBoud == i) {
                m_vTria.push_back(tri);
            }
        }
        m_vBoud[i].iTriaTail = int(m_vTria.size());
    }
}

void TransWorker::fnTetrMaterial()
{
    auto tmpTetr = m_vTetr;
    m_vTetr.clear();

    int iBody = -1, iMedi;
    map<int, int>::iterator iter;
    for (auto &tet : tmpTetr) {
        if (tet.iBody != iBody) {
            iBody = tet.iBody;
            iter = m_mapGeom.find(iBody);
        }
        if (iter == m_mapGeom.end()) continue;
        iMedi = m_vGeom[iter->second].iMedi;
        if (m_vMedi[iMedi].sName == "pec" || m_vMedi[iMedi].sName == "PEC" ||
            m_vMedi[iMedi].dSgm > 1.0e20) {
            iMedi = 0;
        }
        tet.iMedi = iMedi;
        m_vTetr.push_back(tet);
    }
    tmpTetr.clear();
}

void TransWorker::fnLinkTriaTetr()
{
    struct box {
        vector<int> vTet;
    };

    Node r0, r1;
    double dMaxEdge = 0.0, dEdge;
    for (auto &tri : m_vTria) {
        r0 = m_vR[tri.node[1]-1] - m_vR[tri.node[0]-1];
        dEdge = sqrt(r0.x * r0.x + r0.y * r0.y + r0.z * r0.z);
        if (dMaxEdge < dEdge) dMaxEdge = dEdge;
        r0 = m_vR[tri.node[2]-1] - m_vR[tri.node[0]-1];
        dEdge = sqrt(r0.x * r0.x + r0.y * r0.y + r0.z * r0.z);
        if (dMaxEdge < dEdge) dMaxEdge = dEdge;
    }
    dMaxEdge *= 2.0;
    // find the outer box
    r0 = r1 = m_vR[0];
    for (auto &r : m_vR) {
        if (r0.x > r.x) r0.x = r.x;
        if (r0.y > r.y) r0.y = r.y;
        if (r0.z > r.z) r0.z = r.z;
        if (r1.x < r.x) r1.x = r.x;
        if (r1.y < r.y) r1.y = r.y;
        if (r1.z < r.z) r1.z = r.z;
    }
    // construct box
    int iNumX = ceil((r1.x - r0.x) / dMaxEdge);
    int iNumY = ceil((r1.y - r0.y) / dMaxEdge);
    int iNumZ = ceil((r1.z - r0.z) / dMaxEdge);
    map<int, box> boxMap;
    int i, j, k, n;
    for (k = 0; k < iNumZ; ++k) {
        for (j = 0; j < iNumY; ++j) {
            for (i = 0; i < iNumX; ++i) {
                n = (k * iNumX * iNumY + j * iNumX) + i;
                boxMap.insert(pair<int, box>(n, box{}));
            }
        }
    }
    // classify the tet. to each box
    for (size_t m = 0; m < m_vTetr.size(); ++m) {
        r1 = m_vR[m_vTetr[m].node[0]-1];
        for (n = 1; n < 4; ++n) {
            r1 += m_vR[m_vTetr[m].node[n]-1];
        }
        i = ceil((r1.x * 0.25 - r0.x) / dMaxEdge) - 1;
        j = ceil((r1.y * 0.25 - r0.y) / dMaxEdge) - 1;
        k = ceil((r1.z * 0.25 - r0.z) / dMaxEdge) - 1;
        n = (k * iNumX * iNumY + j * iNumX) + i;
        auto iter = boxMap.find(n);
        if (iter != boxMap.end()) iter->second.vTet.push_back(m);
    }
    // check each triangle
    int i1, i2, j1, j2, k1, k2;
    for (auto &tri : m_vTria) {
        r1 = m_vR[tri.node[0]-1];
        for (n = 1; n < 3; ++n) {
            r1 += m_vR[tri.node[n]-1];
        }
        i = ceil((r1.x / 3.0 - r0.x) / dMaxEdge) - 1;
        j = ceil((r1.y / 3.0 - r0.y) / dMaxEdge) - 1;
        k = ceil((r1.z / 3.0 - r0.z) / dMaxEdge) - 1;
        i1 = max(i-1, 0);
        j1 = max(j-1, 0);
        k1 = max(k-1, 0);
        i2 = min(i+1, iNumX-1);
        j2 = min(j+1, iNumY-1);
        k2 = min(k+1, iNumZ-1);
        tri.iTetr = -1;
        for (i = i1; i <= i2; ++i) {
            for (j = j1; j <= j2; ++j) {
                for (k = k1; k <= k2; ++k) {
                    n = (k * iNumX * iNumY + j * iNumX) + i;
                    auto iter = boxMap.find(n);
                    if (iter == boxMap.end()) continue;
                    for (int &m : iter->second.vTet) {
                        n = 0;
                        if (tri.node[0] == m_vTetr[m].node[0] ||
                            tri.node[0] == m_vTetr[m].node[1] ||
                            tri.node[0] == m_vTetr[m].node[2] ||
                            tri.node[0] == m_vTetr[m].node[3]){
                            ++n;
                        }
                        if (tri.node[1] == m_vTetr[m].node[0] ||
                            tri.node[1] == m_vTetr[m].node[1] ||
                            tri.node[1] == m_vTetr[m].node[2] ||
                            tri.node[1] == m_vTetr[m].node[3]){
                            ++n;
                        }
                        if (tri.node[2] == m_vTetr[m].node[0] ||
                            tri.node[2] == m_vTetr[m].node[1] ||
                            tri.node[2] == m_vTetr[m].node[2] ||
                            tri.node[2] == m_vTetr[m].node[3]){
                            ++n;
                        }
                        if (3 == n) {
                            tri.iTetr = m + 1;
                            break;
                        }
                    }
                    if (tri.iTetr > 0) break;
                }
                if (tri.iTetr > 0) break;
            }
            if (tri.iTetr > 0) break;
        }
    }
}

void TransWorker::fnGetWireFrame(PortInfo &port, vector<int> &vWire)
{
    int iHead = 0, iTail = 0;

    vWire.clear();
    for (auto &boud : m_vBoud) {
        if (port.id == boud.id) {
            iHead = boud.iTriaHead;
            iTail = boud.iTriaTail;
            break;
        }
    }
    if (iHead >= iTail) return;

    // get the normal vector of port face
    int i = m_vTria[iHead].node[0] - 1;
    int j = m_vTria[iHead].node[1] - 1;
    int k = m_vTria[iHead].node[2] - 1;
    Node r1, r2, norm;
    r1 = m_vR[j] - m_vR[i];
    r2 = m_vR[k] - m_vR[i];
    cross_mult(r1, r2, port.norm);
    normalize(port.norm);

    // get all edge
    vector<vector<int>> vvEdge;
    auto vEdge = vector<int>(2);
    int iNode;
    for (i = iHead; i < iTail; ++i) {
        for (j = 0; j < 3; ++j) {
            vEdge[0] = m_vTria[i].node[j] - 1;
            vEdge[1] = m_vTria[i].node[(j+1)%3] - 1;
            if (vEdge[0] > vEdge[1]) {
               iNode = vEdge[0];
               vEdge[0] = vEdge[1];
               vEdge[1] = iNode;
            }
            vvEdge.push_back(vEdge);
        }
    }
    sort(vvEdge.begin(), vvEdge.end(), CompareEdge);
    // get all out boundary edge
    auto vvTmp = vvEdge;
    vvEdge.clear();
    k = (int)vvTmp.size() - 1;
    for (i = 0; i < vvTmp.size(); ++i) {
        if (i > 0 && vvTmp[i][0] == vvTmp[i - 1][0] && vvTmp[i][1] == vvTmp[i - 1][1]) continue;
        if (i < k && vvTmp[i][0] == vvTmp[i + 1][0] && vvTmp[i][1] == vvTmp[i + 1][1]) continue;
        vvEdge.push_back(vvTmp[i]);
    }
    vvTmp.clear();
    // extract wire frame
    vEdge.clear();
    vEdge.push_back(vvEdge[0][0]);
    vEdge.push_back(vvEdge[0][1]);
    vvEdge[0].clear();
    iNode = vEdge[1];
    while (true) {
        bool bOut = true;
        for (auto &edge : vvEdge){
            if (edge.empty()) continue;
            if (edge[0] == iNode) {
                iNode = edge[1];
                vEdge.push_back(iNode);
                bOut = false;
                edge.clear();
                break;
            }
            if (edge[1] == iNode) {
                iNode = edge[0];
                vEdge.push_back(iNode);
                bOut = false;
                edge.clear();
                break;
            }
        }
        if (bOut) break;
    }
    if (vEdge[0] == vEdge.back()) vEdge.pop_back();
    // remove aaditional point
    double alph = sin(5 * 3.141592653 / 180.0);
    k = (int)vEdge.size();
    for (i = 0; i < vEdge.size(); ++i) {
        r1 = m_vR[vEdge[(i + 1) % k]] - m_vR[vEdge[i]];
        r2 = m_vR[vEdge[(i - 1 + k) % k]] - m_vR[vEdge[i]];
        normalize(r1);
        normalize(r2);
        cross_mult(r1, r2, norm);
        if (get_length(norm) > alph) vWire.push_back(vEdge[i]);
    }
    // get the center point of port
    r1 = m_vR[vWire[0]];
    for (size_t i = 1; i < vWire.size(); ++i) {
        r1 += m_vR[vWire[i]];
    }
    port.rc = r1 / vWire.size();
    if (abs(port.rc.x) < 1.0e-12) port.rc.x = 0;
    if (abs(port.rc.y) < 1.0e-12) port.rc.y = 0;
    if (abs(port.rc.z) < 1.0e-12) port.rc.z = 0;
    if (abs(port.norm.x) < 1.0e-12) port.norm.x = 0;
    if (abs(port.norm.y) < 1.0e-12) port.norm.y = 0;
    if (abs(port.norm.z) < 1.0e-12) port.norm.z = 0;

    auto tmp = vvEdge;
    vvEdge.clear();
    for (auto &vE : tmp) {
        if (!vE.empty()) vvEdge.push_back(vE);
    }
    if (vvEdge.empty()) return;
    // extract wire frame
    vEdge.clear();
    vEdge.push_back(vvEdge[0][0]);
    vEdge.push_back(vvEdge[0][1]);
    vvEdge[0].clear();
    iNode = vEdge[1];
    while (true) {
        bool bOut = true;
        for (auto &edge : vvEdge){
            if (edge.empty()) continue;
            if (edge[0] == iNode) {
                iNode = edge[1];
                vEdge.push_back(iNode);
                bOut = false;
                edge.clear();
                break;
            }
            if (edge[1] == iNode) {
                iNode = edge[0];
                vEdge.push_back(iNode);
                bOut = false;
                edge.clear();
                break;
            }
        }
        if (bOut) break;
    }
    if (vEdge[0] == vEdge.back()) vEdge.pop_back();
    // remove aaditional point
    k = (int)vEdge.size();
    for (size_t i = 0; i < vEdge.size(); ++i) {
        r1 = m_vR[vEdge[(i + 1) % k]] - m_vR[vEdge[i]];
        r2 = m_vR[vEdge[(i - 1 + k) % k]] - m_vR[vEdge[i]];
        normalize(r1);
        normalize(r2);
        cross_mult(r1, r2, norm);
        if (get_length(norm) > alph) vWire.push_back(vEdge[i]);
    }
}

bool TransWorker::fnOutputMsh(string &sName)
{
    auto indx = sName.find_last_of('.');
    auto sFile = sName.substr(0, indx) + ".msh";
    ofstream fp(sFile);
    if (!fp.is_open()) {
        cout << "Error: Create msh mesh file failed." << endl;
        return false;
    }

    fp << "MESH dimension 3 ElemType Tetrahedra Nnode 4" << endl;
    fp << "Coordinates" << endl;
    int i = 0;
    for (auto &r : m_vR) {
        ++i;
        fp << i << " " << setiosflags(ios::scientific) << setprecision(15) << r.x << " " << r.y << " " << r.z << endl;
    }
    fp << "End Coordinates" << endl << endl;
    fp << "Elements" << endl;
    i = 0;
    for (auto &tet : m_vTetr) {
        ++i;
        fp << i << " " << tet.node[0] << " " << tet.node[1] << " ";
        fp << tet.node[2] << " " << tet.node[3] << endl;
    }
    fp << "End Elements" << endl;
    fp.close();
    return true;
}

bool TransWorker::fnOutputTet(string &sName)
{
    auto indx = sName.find_last_of('.');
    auto sFile = sName.substr(0, indx) + ".tet";
    ofstream fp(sFile);
    if (!fp.is_open()) {
        cout << "Error: Create tet mesh file failed." << endl;
        return false;
    }

    fp << m_vR.size() << " " << m_vTetr.size() << " " << m_vBoud.size() << endl;
    // output node coordinate
    for (auto &r : m_vR) {
        fp << setiosflags(ios::scientific) << setprecision(10) << r.x << " " << r.y << " " << r.z << endl;
    }
    // output element
    for (auto &tet : m_vTetr) {
        fp << tet.node[0] << " " << tet.node[1] << " " << tet.node[2] << " ";
        fp << tet.node[3] << " " << tet.iBody /*tet.iMedi*/ << endl;
    }
    // output boundary
    int iNumTria;
    for (auto &boud : m_vBoud) {
        iNumTria = boud.iTriaTail - boud.iTriaHead;
        fp << boud.sType << " " << boud.id << " " << iNumTria << endl;
        for (int i = boud.iTriaHead; i < boud.iTriaTail; ++i) {
            fp << m_vTria[i].iTetr << " " << m_vTria[i].node[0] << " " ;
            fp << m_vTria[i].node[1] << " " << m_vTria[i].node[2] << endl;
        }
    }
    fp.close();
    return true;
}

bool TransWorker::fnOutputXml(string &sName)
{
    XMLDocument doc;
    const char *head = "<?xml version=\"1.0\" encoding=\"utf-8\"?>";

    doc.Parse(head);
    auto pRoot = doc.NewElement("CEMS");
    auto pFreq = doc.NewElement("Frequency");
    pFreq->SetAttribute("Unit", m_para.sFreqUnit.c_str());
    pFreq->SetAttribute("AnalysisFreq", m_para.dFreq0);
    if (m_para.iNumFreq > 0) {
        pFreq->SetAttribute("MinFreq", m_para.dMinFreq);
        pFreq->SetAttribute("MaxFreq", m_para.dMaxFreq);
        pFreq->SetAttribute("NumFreq", m_para.iNumFreq);
        pFreq->SetAttribute("SweepType", m_para.sSwpType.c_str()); // Discrete, Interpolate, Fast
    }
    pRoot->InsertEndChild(pFreq);

    auto pSymm = doc.NewElement("Symmetry");
    pSymm->SetAttribute("YOZ", 0);
    pSymm->SetAttribute("XOZ", 0);
    pSymm->SetAttribute("XOY", 0);
    pRoot->InsertEndChild(pSymm);

    auto pMedi = doc.NewElement("Material");
    pMedi->SetAttribute("Count", static_cast<int>(m_vMedi.size()-1));
    pMedi->SetAttribute("FreqRevise", "false");
    for (unsigned i = 1; i < m_vMedi.size(); ++i) {
        auto pSon = doc.NewElement(m_vMedi[i].sName.c_str());
        pSon->SetAttribute("ID", i);
        if (m_vMedi[i].iTensor == 0) {
            pSon->SetAttribute("Tensor", "false");
            pSon->SetAttribute("Epsilon", m_vMedi[i].dEps);
            pSon->SetAttribute("tanE", m_vMedi[i].dTan);
            pSon->SetAttribute("Mu", m_vMedi[i].dMiu);
            pSon->SetAttribute("tanM", m_vMedi[i].dTanM);
            pSon->SetAttribute("Sigma", m_vMedi[i].dSgm);
        } else {
            pSon->SetAttribute("Tensor", "true");
            auto pRealEps = doc.NewElement("RealEps");
            for (int idx1 = 0; idx1 < 3; ++idx1) {
                for (int idx2 = 0; idx2 < 3; ++idx2) {
                    pRealEps->SetAttribute(("e" + to_string(idx1 + 1) + to_string(idx2 + 1)).c_str(), m_vMedi[i].vEpsRe[idx1][idx2]);
                }
            }
            pSon->InsertEndChild(pRealEps);
            auto pImagEps = doc.NewElement("ImagEps");
            for (int idx1 = 0; idx1 < 3; ++idx1) {
                for (int idx2 = 0; idx2 < 3; ++idx2) {
                    pImagEps->SetAttribute(("e" + to_string(idx1 + 1) + to_string(idx2 + 1)).c_str(), m_vMedi[i].vEpsIm[idx1][idx2]);
                }
            }
            pSon->InsertEndChild(pImagEps);
            auto pRealMiu = doc.NewElement("RealMiu");
            for (int idx1 = 0; idx1 < 3; ++idx1) {
                for (int idx2 = 0; idx2 < 3; ++idx2) {
                    pRealMiu->SetAttribute(("u" + to_string(idx1 + 1) + to_string(idx2 + 1)).c_str(), m_vMedi[i].vMiuRe[idx1][idx2]);
                }
            }
            pSon->InsertEndChild(pRealMiu);
            auto pImagMiu = doc.NewElement("ImagMiu");
            for (int idx1 = 0; idx1 < 3; ++idx1) {
                for (int idx2 = 0; idx2 < 3; ++idx2) {
                    pImagMiu->SetAttribute(("u" + to_string(idx1 + 1) + to_string(idx2 + 1)).c_str(), m_vMedi[i].vMiuIm[idx1][idx2]);
                }
            }
            pSon->InsertEndChild(pImagMiu);
        }
        pMedi->InsertEndChild(pSon);
    }
    pRoot->InsertEndChild(pMedi);

    auto pPost = doc.NewElement("PostSet");
    auto pFar = doc.NewElement("FarField");
    pFar->SetAttribute("Count", 1);
    auto pSon = doc.NewElement("sphere");
    pSon->SetAttribute("NumPhi", 181);
    pSon->SetAttribute("MinPhi", 0);
    pSon->SetAttribute("MaxPhi", 180);
    pSon->SetAttribute("NumTheta", 361);
    pSon->SetAttribute("MinTheta", -180);
    pSon->SetAttribute("MaxTheta", 180);
    pFar->InsertEndChild(pSon);
    pPost->InsertEndChild(pFar);

    auto pNear = doc.NewElement("NearField");
    pNear->SetAttribute("Setup", "false");
    pNear->SetAttribute("Count", 0);
    pPost->InsertEndChild(pNear);
    // output current setting
    auto pCurr = doc.NewElement("Current");
    pCurr->SetAttribute("Setup", "false");
    pPost->InsertEndChild(pCurr);
    pRoot->InsertEndChild(pPost);

    // output the voltage
    auto pVolt = doc.NewElement("Voltage");
    pVolt->SetAttribute("Count", int(m_vLump.size()+m_vPort.size()));
    for (auto &src : m_vLump) {
        auto pSrc = doc.NewElement(src.sName.c_str());
        pSrc->SetAttribute("Group", 1);
        auto pVal = doc.NewElement("Value");
        pVal->SetAttribute("Mag", src.dMag);
        pVal->SetAttribute("Pha", src.dPha);
        pSrc->InsertEndChild(pVal);
        pVolt->InsertEndChild(pSrc);
    }
    for (auto &src : m_vPort) {
        auto pSrc = doc.NewElement(src.sName.c_str());
        pSrc->SetAttribute("Group", 1);
        auto pVal = doc.NewElement("Value");
        pVal->SetAttribute("Mag", src.dMag);
        pVal->SetAttribute("Pha", src.dPha);
        pSrc->InsertEndChild(pVal);
        pVolt->InsertEndChild(pSrc);
    }
    pRoot->InsertEndChild(pVolt);

    auto pSolvr = doc.NewElement("Solver");
    pSolvr->SetAttribute("Parallel", "true");
    pSolvr->SetAttribute("Mpiexe", "C:/program file/MPICH2/bin/mpiexec");
    pSolvr->SetAttribute("NumCore", 4);
    pSolvr->SetAttribute("Integration", "Normal");    // Fine, Normal, Coarse
    pRoot->InsertEndChild(pSolvr);

    auto pDoma = doc.NewElement("Domain");
    pDoma->SetAttribute("Count", 1);
    auto pMod = doc.NewElement("Model");
    pMod->SetAttribute("ModelName", "model");

    int m = 0;
    string sMesh;
    auto indx = sName.find_last_of('\\');
    if (indx != string::npos) m = indx + 1;
    indx = sName.find_last_of('/');
    if (indx != string::npos && indx >= m) m = indx + 1;
    indx = sName.find_last_of('.');
    if (indx != string::npos) {
        sMesh = sName.substr(m, indx - m);
    } else {
        sMesh = sName.substr(m);
    }

    pMod->SetAttribute("MeshFile", sMesh.c_str());
    pMod->SetAttribute("Method", "FEM");
    auto pGeom = doc.NewElement("Geometry");
    pGeom->SetAttribute("Count", static_cast<int>(m_vGeom.size()));
    for (int i = 0; i < m_vGeom.size(); ++i) {
        auto pSon = doc.NewElement("Geom");
        pSon->SetAttribute("Name",m_vGeom[i].sName.c_str());
        pSon->SetAttribute("ID", m_vGeom[i].iObjID);
        pSon->SetAttribute("MaterialID", m_vGeom[i].iMedi);
        pGeom->InsertEndChild(pSon);
    }
    pMod->InsertEndChild(pGeom);
    auto pFEM = doc.NewElement("SetFEM");
    pFEM->SetAttribute("Order", m_para.iOrder);
    pFEM->SetAttribute("SolverType", "Direct"); // Direct, Iteration, DDM
    pFEM->SetAttribute("MaxIterNum", 20);
    pFEM->SetAttribute("MaxIterErr", 0.005);
    // parameter for eigen mode
    pFEM->SetAttribute("SolutionType", "Driven");
    // parameter for adaptive mesh
    if (m_para.iAdapMaxNum > 0) {
        pFEM->SetAttribute("AdaptMinStep", m_para.iAdapMinNum);
        pFEM->SetAttribute("AdaptMaxStep", m_para.iAdapMaxNum);
        pFEM->SetAttribute("AdaptMinPass", m_para.iAdapPassed);
        pFEM->SetAttribute("AdaptDeltaSp", m_para.dAdaptDelta);
        pFEM->SetAttribute("AdaptMeshRatio", m_para.dAdaptRatio);
    }
    pFEM->SetAttribute("OutBoundary", "Radiation");
    if (m_para.bSolvInside) {
        pFEM->SetAttribute("MetalInside", "true");
    }else {
        pFEM->SetAttribute("MetalInside", "false");
    }
    pMod->InsertEndChild(pFEM);

    if (!m_vLump.empty()) {
        auto pPort = doc.NewElement("Lumpport");
        pPort->SetAttribute("Count", static_cast<int>(m_vLump.size()));
        for (auto &port : m_vLump) {
            auto pSon = doc.NewElement(port.sName.c_str());
            pSon->SetAttribute("Type", port.sType.c_str());
            pSon->SetAttribute("BoundID", port.id);
            pSon->SetAttribute("a",  port.a);
            pSon->SetAttribute("b", port.b);
            pSon->SetAttribute("X", port.rc.x);
            pSon->SetAttribute("Y", port.rc.y);
            pSon->SetAttribute("Z", port.rc.z);
            pSon->SetAttribute("Ex", port.ei.x);
            pSon->SetAttribute("Ey", port.ei.y);
            pSon->SetAttribute("Ez", port.ei.z);
            pSon->SetAttribute("Tx", port.norm.x);
            pSon->SetAttribute("Ty", port.norm.y);
            pSon->SetAttribute("Tz", port.norm.z);
            pSon->SetAttribute("Z0", port.dZ0);
            pPort->InsertEndChild(pSon);
        }
        pMod->InsertEndChild(pPort);
    }
    if (!m_vPort.empty()) {
        auto pPort = doc.NewElement("Waveport");
        pPort->SetAttribute("Count", static_cast<int>(m_vPort.size()));
        for (auto &port : m_vPort) {
            auto pSon = doc.NewElement(port.sName.c_str());
            pSon->SetAttribute("Type", port.sType.c_str());
            pSon->SetAttribute("BoundID", port.id);
            pSon->SetAttribute("a",  port.a);
            pSon->SetAttribute("b", port.b);
            pSon->SetAttribute("X", port.rc.x);
            pSon->SetAttribute("Y", port.rc.y);
            pSon->SetAttribute("Z", port.rc.z);
            pSon->SetAttribute("Ex", port.ei.x);
            pSon->SetAttribute("Ey", port.ei.y);
            pSon->SetAttribute("Ez", port.ei.z);
            pSon->SetAttribute("Tx", port.norm.x);
            pSon->SetAttribute("Ty", port.norm.y);
            pSon->SetAttribute("Tz", port.norm.z);
            pSon->SetAttribute("Z0", port.dZ0);
            pPort->InsertEndChild(pSon);
        }
        pMod->InsertEndChild(pPort);
    }
    pDoma->InsertEndChild(pMod);
    pRoot->InsertEndChild(pDoma);

    doc.InsertEndChild(pRoot);

    // get the absolute file name
    indx = sName.find_last_of('.');
    auto sFile = sName.substr(0, indx) + ".xml";
    if (doc.SaveFile(sFile.c_str()) == XML_SUCCESS) {
        return true;
    } else {
        cout << "Error: Create project file failed." << endl;
        return false;
    }
}
